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Abstract 


We  analyze  sequential  detection  for  diffusion  type  signals  both  in  the  fixed  prob¬ 
ability  of  error  formulation  and  in  the  Bayes  Formulation.  We  show  that  the  optimal 
strategy  in  both  cases  is  of  the  threshold  type  with  explicitly  computable  thresholds. 
We  provide  efficient  numerical  schemes  for  computing  approximations  to  the  likeli¬ 
hood  ratio  and  provide  an  implementation  via  a  special  purpose  VLSI  processor  for 
real  time  processing  in  the  scalar  diffusion  case.  Finally,  we  describe  the  DELPHI 
expert  system  which  is  under  development.  Its  purpose  is  to  provide  an  integrated 
system  level  design  tool  for  sequential  real  time  detection  and  estimation. 


Introduction 


In  the  present  paper  we  analyze  in  detail  the  sequential  detection  problem  for  dif¬ 
fusion  process  signals.  In  particular  we  address  questions  related  to  optimal  strategies, 
numerical  computation  of  the  sufficient  statistics,  real-time  architectures  for  imple¬ 
menting  the  resulting  algorithms,  and  describes  a  sophisticated  system  level  design 
tool:  the  DELPHI  system.  The  techniques  presented  here  can  be  extended  to  apply 
to  more  general  signal,  observation  models  such  as  multi-dimensional  diffusions,  jump 
point  process  models  and  mixed  models.  We  shall  not  address  them  here. 

While  Gaussian  signal  processing  theory  and  design  have  reached  a  certain  de¬ 
gree  of  completeness,  the  corresponding  developments  for  non-Gaussian  signal  pro¬ 
cessing  are  unsatisfactory  from  the  designer’s  point  of  view.  The  main  reason  is  that 
no  systematic  effort  has  been  undertaken  to  transform  the  theoretical  advances  in 
non-Gaussian  detection  and  estimation  theory  into  practical  design  methods.  In  the 
meantime  the  engineering  specifications  on  signal  processors  are  becoming  increas¬ 
ingly  more  demanding,  resulting  in  an  unavoidable  ascending  degree  of  complexity. 

Recent  progress  in  stochastic  processes,  and  in  particular  in  non-Gaussian  detec¬ 
tion  and  nonlinear  altering  theory,  hold  great  promise.  It  is  well  known  that  every 
signal  detection  and  classification  problem  requires  the  solution  of  some  underlying 
estimation  problem  (Kailath  [1969]).  Despite  the  fact  that  this  has  been  known 
for  many  years,  its  impact  on  the  resolution  of  complex  signal  processing  problems 
in  non-Gaussian  environments  has  been  very  limited  (see  the  references  of  (Kailath 
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[1969])  for  some  ingenious  applications).  The  primary  reason  for  the  limited  appli¬ 
cability  of  the  powerful  machinery  of  nonlinear  estimation  in  signal  detection  and 
classification  problems  has  been  the  formidable  computational  complexity  called  for 
by  these  theories:  at  least  a  stochastic  nonlinear  partial  differential  equation  has 
to  be  solved  on-line.  However,  recently  three  events  have  caused  a  serious  reexam¬ 
ination  of  these  issues.  First,  recent  progress  in  nonlinear  estimation  succeeded  in 
replacing  the  fundamental  stochastic  nonlinear  partial  differential  equation  with  a 
linear  non-stochastic  partial  differential  equation,  (Davis  [1981])  the  so  called  robust 
version  of  the  Zakai  equation.  Second,  the  advent  of  systolic  VLSI  arrays  (as  well 
as  other  special  purpose  VLSI  array  processors)  led  to  the  proposal  of  a  design  of  a 
VLSI  processor  that  could  implement  the  Zakai  equation  in  real  time  (Baras  [1981]). 
Third,  the  maturation  of  artihcial  intelligence  has  provided  a  plethora  of  design  tools 
whose  impact  on  other  disciplines  has  not  even  been  tested  yet. 

The  nonlinear  filtering  problem,  in  its  various  manifestations  (depending  on  sig¬ 
nal  and  observation  models),  is  at  the  heart  of  many  interesting  and  significant  signal 
processing  problems.  We  have  identified,  among  others,  digital  signal  processing 
(such  as  pulse  amplitude  modulation,  delta  modulation,  adaptive  delta  modulation, 
and  speech  processing),  direction  finding  receivers,  digital  phase  lock  loops,  adaptive 
sonar  and  radar  arrays,  sequential  detection,  simultaneous  detection  and  estimation, 
sensor  scheduling  in  data  fusion  problems,  adaptive  stochastic  control,  stochastic  con¬ 
trol  with  partial  observations,  and  nonlinear  observers. 

We  next  offer  some  justification  for  the  selection  of  the  research  problem  and 
technical  approach.  First  regarding  VLSI  architectures,  we  note  that  recently  a  lot  of 
attention  and  research  effort  has  been  devoted  to  VLSI  architectures  for  linear  signal 
processing  schemes  (Kailath  Whitehead  ]).  However,  the  need  for  special  architectures 
is  more  dramatic  in  the  case  of  nonlinear  signal  processors,  such  as  needed  in  non- 
Gaussian  problems.  To  focus,  we  recall  that  the  sequential  detection  problem  or  the 
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nonlinear  filtering  problem,  calls  for  a  processor  that  operates  on  data  in  “real  time” . 
What  is  meant  by  “real  time”?  Since  it  is  clear  that  we  are  led  to  some  sort  of 
digital  circuit  implementation,  real  time  means  that  the  allowable  sampling  period 
for  the  observed  data  must  be  larger  that  the  computational  cycle  for  the  digital 
device  implementing  the  algorithm.  In  the  case  of  the  nonlinear  filtering  problem, 
the  theoretical  analysis  produces  an  algorithm  that  calls  for  the  solution  of  a  stochastic 
partial  differential  equation.  Hence  a  real  time  processor  based  on  this  algorithm  must 
provide  fast  approximations  to  the  solution  of  this  equation.  To  a  great  extent  our 
work  has  been  initiated  by  the  desire  for  developing  a  theory  and  design  to  account 
for  real  time  implementation  constraints  for  non-Gaussian  detection  and  non-linear 
filtering. 

The  structure  of  this  paper  is  as  follows. 

In  Section  1,  we  provide  a  detailed  analysis  of  the  sequential  detection  problem 
for  diffusion  signal  and  observation  models.  Both  the  Bayesian  and  fixed  probability 
of  error  formulations  are  examined.  We  show  that  in  both  cases  the  optimal  strategy 
is  of  threshold  type  with  thresholds  that  can  be  computed  explicitly. 

In  Section  2,  we  provide  a  detailed  treatment  of  the  numerical  analysis  of  the 
Zakai  equation  using  semigroup  methods.  A  general  approximation  theorem  is  pre¬ 
sented.  This  theorem  is  intended  to  be  used  to  check  convergence  of  individual  ap¬ 
proximation  schemes. 

In  Section  3,  we  describe  a  VLSI  architecture  that  can  implement  the  algorithms 
presented  in  Section  2  in  real  time.  Initial  estimates  indicate  that  for  100  mesh  points, 
we  can  perform  5000  calculations  per  second. 

In  Section  4,  we  describe  the  DELPHI  expert  system,  which  is  under  development 
at  Maryland.  It  affords  an  engineer  a  sophisticated  and  “  user  friendly”  design  tool 
for  real-time  circuit  design  iin  non-Gaussian  detection  and  estimation  problems. 

Finally,  we  discuss  briefly  some  future  direction  of  the  work  described  here. 


1.  Sequential  Detection  of  Diffusion  Signals 


1.1  Introduction 

In  this  chapter,  we  consider  the  binary  sequential  hypothesis  testing  problem.  Here, 
we  are  given  an  JR”- valued  signal  process  {xt,  t  >  0}  which  satisfies  the  stochastic 
differential  equation 


dxt  =  f{xt)  dt  +  g{xt)  dwf 

Xo  =  ^ 


(1.1) 


where  {wt,  t  >  0}  is  an  R'^-valued  standard  Brownian  motion.  However,  we  cannot 
observe  {xt,  t  >  0}  directly,  instead  we  only  observe  the  increments  dyt  of  an  Re¬ 
valued  stochastic  process  {yt,  t  >  0}.  Under  each  hypothesis  the  observed  data  is 
the  output  of  a  stochastic  differential  equation,  i.e.. 

Under  Hi  :  dyt  =  h{xt)  dt  +  dvt 


(1.2) 


Under  Hq  :  dyt  =  dvt 

where  {vt,  t  >  0}  is  an  R^-valued  standard  Brownian  motion  which  is  independent 
of  {tvt,  t  >  0}. 


We  assume  that  for  all  x,  y  in  R”,  the  functions  /,  g,  and  h  satisfy  the  Lipschitz 
condition. 


ll/(a:)  -  f{y]\\  +  ||ff{a:)  -  ?(y)l|  +  i|/i(a:)  -  h{y)\\  <  K\\x  -  y\\  (1.3) 


and  the  growth  condition 

ii/wii" + ikwr + iiAwii"  <  + ikin-  (1-4) 
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These  conditions  guarantee  that  the  stochastic  differential  equations  in  (l.l)  have 
unique  continuous  strong  solutions  (Arnold  [1974]). 

Data  is  observed  continuously  starting  at  an  initial  time  which  is  taken  for  con¬ 
venience  to  be  zero.  At  each  time  f  >  0,  the  decision-maker  can  either  declare  one  of 
the  hypotheses  to  be  true  or  continue  collecting  data.  The  decision-maker  selects  his 
decision  based  on  the  data  collected  up  to  time  t,  so  as  to  minimize  an  appropriate 
cost  function. 

We  will  present  both  the  fixed  probability  of  error  and  the  Bayesian  formulations 
for  this  problem.  In  both  formulations,  we  have  a  measurable  space  (fl,  J),  on  which 
we  are  given  two  probability  measures  Po,  Pi,  and  the  random  process  {yt,  t  >  0}. 
When  hypothesis  Hq  (respectively  Hi)  is  valid  the  statistics  of  the  observed  process 
{Vt,  ^  ^  0}  are  governed  by  measure  Pq  (respectively  Pi). 

The  difference  between  these  formulations  is  how  they  prescribe  the  cost  function 
to  be  minimized.  For  each  formulation,  we  show  that  the  appropriate  cost  function  is 
minimized  by  a  threshold  policy.  The  essential  difference  between  these  formulations 
is  how  they  prescribe  the  thresholds. 

More  precisely,  a  decision  policy  involves  the  selection  of  a  termination  time  r, 
and  of  a  binary  valued  decision  6.  If  ^  =  1,  we  shall  accept  hypothesis  Pi;  if  ^  =  0  we 
shall  accept  hypothesis  Hq.  Let  //  denote  the  cr-algebra  generated  hy  {y 3,  s  <t}. 

Definition  1.1.1.  An  admissible  decision  policy  is  any  pair  u  =  (r,  6)  of  RV’s  where 
r  is  an  J^^'-stopping  time,  and  6  is  an  /,?'-measurable  {0,  l}-valued  RV.  The  collection 
of  all  admissible  decision  policies  will  be  denoted  by  U. 

Definition  1.1.2.  A  policy  u  in  1/  is  a  threshold  policy  or  of  threshold  type  if  there 
exists  constants  A  and  B,  with  O<A<I<P<00  and  A  ^  B,  such  that 


T  =  inf(t  >  0  I  At  ^  (A,  P)) 


r  1,  A.r>  B 
(  0,  Ar  <  A 


(1.5) 

(1.6) 
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Here  A*  is  the  likelihood  ratio  associated  with  this  problem,  namely 

At  =  exp{j  hjdys-^j  (1.7) 

where  ^  denoted  transpose,  and 

ht  =  Eiih{x)\rn-  (1-8) 

This  chapter  is  organized  as  follows.  First,  we  present  general  results  on  threshold 
policies.  Then,  we  show  how  threshold  policies  solve  the  fixed  probability  of  error 
problem.  Then,  we  show  that,  with  the  cost  function  given  in  (1.42),  threshold 
policies  solve  the  Bayesian  problem.  Finally,  we  state  the  generalization  to  the  case 
where  Hq,  like  Hi,  includes  a  function  of  a  signal  process  {x°,  t  >  O}.  In  both  of 
these  problems,  we  show  how  to  find  the  optimal  threshold  policy. 

Throughout  this  chapter,  we  make  the  following  technical  assumptions. 

(Tl)  £;4lA(a:t)|)  <  oo,  t  >  0. 

(T2)  Pi(;”||fc.l|**  =  oc)  =  l. 

(T3)  <  00,  0  <  f  <  oo. 

where  ht  is  defined  in  (1.8). 
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1.2  Threshold  Policies 

For  each  umU,  let  a[u)  and  /3{u)  be  the  false  alarm  and  miss  probabilities  associated 
with  u,  respectively,  i.e. 


a(tt)  =  Po(^  =  1)  /?(■**)  =  ^i(  <5  =  0).  (1.9) 

Throughout  this  chapter  o;  and  will  be  positive  constants  that  satisfy  the  inequality 

a  +  <  1. 

A  threshold  policy  u  in  U  will  be  described  in  the  form  (l.5)-(l.8)  and  will  be 
identified  with  the  threshold  constants  (A,  B).  Let  F  be  the  collection  of  all  threshold 
policies  in  U. 

We  will  now  prove  some  results  about  threshold  policies.  These  proofs  are  taken 
from  (Liptser  &;  Shiryayev  [1978,  Chapter  17,6]). 

Lemma  1.2.1.  For  a  threshold  policy  u  in  F, 

Pi(T  <oo)  =  l,  i  =  0,l.  (1.10) 


Proof:  We  will  show  the  result  for  i  =  1  as  a  similar  argument  works  for  the  case 
i  =  0.  Let  On  be  the  sequence  of  //-stopping  times  defined  by 

o„  =  inf(t>0|  /  ||^s|/ds>n) 

Jo 

for  each  n  =  0, 1, . . ..  The  very  definition  of  r  yields  the  inequalities 

rrhcTn  ^  j  prAan  ^ 

A  <  exp(  /  hj  dya  —  -  Insll^  ^ 

Jo  2  Jq 


or  equivalently 


log  A 


rTAcr„ 


hj  dys 


IL 


ds  <  log  B. 


(1.11) 
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Under  Pi,  it  can  be  shown  (Liptser  Shiryayev  [1977])  that  {yt,  t  >  0}  satisfies  the 
stochastic  differential  equation 

dyt  =  htdt  +  dW  t 


with  {W t,  t  >  0}  a  standard  Brownian  motion,  and  upon  substituting  (1.12)  into 
(l.ll)  we  see  that 


logA<  /  hJdWa  +  -  /  ||fes||^  ds  <  log  jB  (1-13) 

Jo  2  Jq 

From  the  definition  of  On,  Jq”'  H^sjP  ds  —  n  ,  from  which  it  follows  that 

I'-rAan  _ 

Eii  hJdWa)=0 

Jo 

Now,  by  taking  the  expectation  of  (1.13)  under  Pi,  we  get  from  (1.14)  that 

1  frAcTn  ^  r-rAcrn  ^  _ 

Pl(logA,A<rJ=Pl(i  /  \\ha\\^ds+  /  h^W s) 

2  Jo  Jo 

1  ^ 

=  Ei{-  ||;i3ird5)<logP. 

^  Jo 


(1.14) 


Since  (1.14)  is  true  for  all  n  and  <t„  f  oo  it  follows  that 

Eiil  \\k\?  ds) 


<  oo 


(1.15) 


whence 


oo  > 


ilfc.ir*}. 


The  conclusion 


Pi  (r  =  oo)  =  0 


is  now  readily  obtained  via  (T2). 

Since  At  has  a.s.  continuous  sample  paths  and  since  r  is  a.s.  finite,  we  can  conclude 
that  for  threshold  policies,  A^  takes  on  the  values  A  or  B  (Pq-  and  Pi -a.s.). 
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Lemma  1.2.2.  For  a  threshold  policy  u  in  F  with  thresholds  A  and  B, 

1-A  .  A[B-1) 


a{u)  — 


B-A' 


13  (u) 


B-A 


(1.16) 


Proof:  From  the  comment  following  Lemma  1.2.1,  we  know  that 

o:(tt)=Po(^=l)  =  Po(A^  =  B)  (1.17) 

and 

I3{u)  =  PiiS  -  0)  =  Pi  (A,  =  A)  (1.18) 

From  Girsanov’s  Theorem  (Liptser  &  Shiryayev  [1977])  we  see  that 

hence 

Pi(A,  =  A)  =  Po(1(a.=a)A.)  =  APo(A,  =  A)  (1.19) 

similarly, 

Po(A,  =  P)  =  Pi(  1(a^=b)A7'  )  =  ^Pi(A.  =  B).  (1.20) 

It  also  follows  from  the  comment  following  Lemma  1.2.1  that 

Po(A,  =  A)  =  1-Po(A,  =  P)  (1.21) 

and 

Pi(A^  =  P)  =  l-Pi(A,  =  A).  (1.22) 

From  (1.19)-(1.22)  and  some  simple  algebra,  we  get  (1.16)  □ 

This  result  has  several  immediate  consequences. 

Corollary  1.2.3.  For  any  threshold  policy  u  in  F, 


«(«)  +  ^(u)  <  1. 


10 


Proof:  This  follows  from  (1.16)  above  and  from  the  fact  that  O<A<I<P<00 
with  A  ^  B. 


Corollary  1.2.4.  Let  u  be  a  threshold  policy  in  F  with  A  and  B  defined  by 


A  = 


1  —  Q! 


B  = 


1-/3 


OL 


(1.23) 


where  a  +  /3  <  1,  then 


Q!(li)  =  Q!  ^{u)  =  /?. 


Proof:  This  is  obtained  by  direct  substitution  of  (1.23)  into  (I.I6). 

Lemma  1.2.5.  Let  u  =  (r,  ^)  be  any  policy  in  U  with  oi{u)  +  ^[u]  <  1  and  let 


a  :=  Q!(u)  /?  :=  /3{u). 


Define  u*  =  {t*,6*)  to  be  the  threshold  policy  in  F  with  parameters  (A*,B*)  that 
correspond  to  the  pair  (o,/3)  as  defined  in  (1.23). 

Then 


Eq 

El 


llb^ll^ds)  =  2w{a,/3) 
llb^ll^ds)  =  2w{P,a) 


(1.24) 

(1.25) 


where 


«;(a:,y)  ;=  (1  -  x)  log - l-a^log- - ,  0  <  <  1. 

y  1-  y 


(1.26) 


Proof:  We  will  only  show  (1.24)  since  a  similar  argument  works  for  (1.25).  Let 

Lt  =  log  At  =  f  hj  dys  -  ~  f  1|L||^  ds,  t  >  0. 

Jo  ^  Jo 

Consider  the  solution  gi{x)  of  the  boundary  valued  problem 

.  i  =  0, 1  (1.27) 

ffi(log  A)  =  E)  =  0. 
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on  the  interval  [log  A,  log  B] .  Elementary  calculations  show  that 

(R  _  4  R#>-*  R  \ 

B  A  +  (1-28) 

go(a:)  =  2  ^  log  ^  -  log  B  ,  (1.29) 

for  log  A  <  X  <  logjB.  From  (1.15),  (1.26),  (1.28),  and  (1.29),  we  now  see  that 

go{0)  =  2w{a,^)  (1.30) 

gi(0)  =  2«;(;d,Q;).  (1.31) 

By  applying  Ito’s  rule  to  go(.&TA<T„),  we  get 

/•rA<T„  ^  j  pTAcTn  .  ,  ^ 

9o{^TA(Tn)  ~  50(0)  “h  J  SoiI^s)hg  dyg  d”  2  y  (■^«)  ~  9o{Ls)j  |l^s||  ds 

prACn  ^  prAcTn  ^ 

=  ?o(0)+  /  g'o{Lg)h'^  dyg  -  (1-32) 

Jo  Jo 

Again  from  the  definition  of  Cn  and  the  fact  that  go(x)  is  bounded  for  x  6 
[logA,  logB],  we  see  that 

rTAUn  ^  rTAcTn  ^ 

Eo{  /  gQ{Ls]hJ  dyg  )^Eo{  9o{Ls)hJ  dvg)=0 

Jo  Jo 

and  consequently  upon  taking  the  expectation  of  (1.32)  with  respect  to  Pq  we  get 

nTAcTn  ^ 

Eo{9o{LrAffn))  —  9o{^)  —  Eo{  jj^all^ds).  (1.33) 

Jo 

Letting  n  go  to  infinity  in  (1.33),  we  finally  obtain 

0  =  go{0)-Eo{f  \\hg\\^ds), 

Jo 

whence  (1.24).  The  equality  (1.25)  can  be  established  using  similar  methods.  □ 


12 


1.3  Fixed  Probability  of  Error  Formulation 

Given  0  <  a, /?  <  1  with  o  +  /3  <  1,  let  U{a,^)  be  the  set  of  all  admissible  policies  u 
in  U  such  that 

Q'(tt)  <  O'  j3(u)  <  13 


The  fixed  probability  of  error  formulation  to  the  sequential  hypothesis  testing 
problem  requires  the  solution  of  the  following. 

Problem  (Pf)-  Find  u*  in  U(a,/3)  such  that  for  all  u  in 

U{a,P), 


Ei{  r  \\hs\\^  ds)  >  Ei{  r  ll^.ll^da),  i  = 
Jo  Jq 


(1.34) 


Theorem  1.3.1.  If  u*  is  the  threshold  policy  with  constants  {A*,B*)  defined  by 


then  u*  solves  problem  (Pp )• 

Proof;  From  Lemma  1.2.2,  it  follows  that  u  G  ll{a,l3).  Hence,  we  only  need  to  show 
u*  is  optimal  is  the  sense  of  (1.33).  To  this  end,  let  u  =  (r,  S)  be  any  policy  in  U  (a,  /3). 
From  Lemma  1.2.1  we  see  that 

II^.IP  *)  =  Ei(£  a?  dv.  -  i  £  Ijfc.ll"  da) 

=  £;i(logAt) 

=  -EiilogA-^) 

and  from  Jensen’s  inequality 

Eiil  [  ^  -£’i(log  A7I)  (1.35) 

=  -P^{8  =  l).Ei(logA7^  I  5  =  1)  -  Fi(^  =  0)fi;i(logA7^  I  (5  ==  0) 

>  -Fi(^  =  l)log£;i(A7^  I  5  =  1)  -Pi(^  =  0)logEi(A;'  K  =  0). 
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Since 


Pq(^S  =  i)  —  Ei{l(^s=i)^T 

=  Ei{l^s=i)E,{K"\S  =  i)) 

=  Pii8  =  i)Ei{K"\8  =  i),  (1.36) 


the  inequality  (1.35)  becomes 


El 


=  (1  -  /3)  log  ^  +  /?  log  ^ 


O! 


1  —  O' 


The  inequality 


ds). 


(1.37) 


Eo{r\\kfds)>Eo{r  iih.ii^ds) 

Jo  Jo 


can  be  established  in  a  similar  fashion. 
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1.4  Bayesian  Formulation 

For  the  Bayesian  formulation,  let  H  be  an  {0,  l}-valued  RV  indicating  the  true  hy¬ 
pothesis.  By  (p  we  denote  the  a  priori  probability  that  hypothesis  Hi  is  true.  We 
consider  a  probability  measure  P  on  (O,  J)  such  that 

P{H=l)  =  >p  P{H  =  0)^l-ip  (1.38) 

and  such  that  for  every  A  G  I 

P{A)  =  <pPi{A)  +  {l-p)PQiA),  (1.39) 


where  Pi  and  Pq  are  the  measures  defined  in  Section  1.1. 

We  shall  assume  the  cost  of  observation  to  accrue  according  to  k  ||hs|pds, 
where  k  >  0  and  {ht,  t  >  0}  is  defined  by  (1.8).  The  average  cost  due  to  data 
collection  is  simply 


ji{T)=E{k  r  whtw^dt) 

Jo 


(1.40) 


where  the  expectation  is  taken  under  P.  The  costs  associated  with  the  binary  decision 


6  are  given  by 


C{H,6) 


Cl,  when  H  =  1  and  ^  =  0; 
C2,  when  H  —  0  and  <5=1; 
0,  otherwise. 


where  ci  >  0  and  C2  >  0.  The  average  cost  due  to  the  selection  of  6  is 


(1.41) 


J2{6)=E[C{H,6)] 

=  ciP{H  =  1,6  =  0)  +  C2P[H  =  0, «  =  1) 

=  CiipPi{S  =  0)  +  C2(l  -  ‘p)Po{6  =  1).  (1-42) 

If  u  =  (’■)  ^)  is  any  admissible  policy,  then  the  corresponding  average  expected  cost  is 

J(u)  =  E(k  r  l\h,jj^ds  +  C(H,S)).  (1.43) 

Jo 

The  Bayesian  approach  to  sequential  detection  requires  the  solution  of  the  fol¬ 


lowing  optimization  problem. 
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Problem  (Pb):  Given  ip  G  (0, 1),  find  u*  in  U  such  that, 


J{u*)=  inf  J(u). 

u€U 


(1.44) 


Let  TTt  =  P[H  =  1  I  be  the  a  posteriori  probability  of  the  hypothesis  Hi 
given  7^.  From  the  definition  of  irt  we  see  that 

P{H=l\in  =  E{H\in 


=  pE{ 


dPi 


p- 


d{pPi  +  (1  -  ^)Po) 


d{pPi  +  (1  -  <^)Po) 

From  Girsanov’s  Theorem,  we  know  that 

dPi 


(1.45) 


dPo 


(//)=At,  t>0, 


hence  (1.45)  becomes 


^t  =  P(^=l|^/) 


pAt 


pAt  +  (1  -  p)' 

Lemma  1.4.1.  Let  u  =  (r,  ^)  be  any  policy  in  U  with  a{u)  +  /3(u)  <  1  and  pose 

a  :=  a{u)  /3  :=  /3(u). 

If  u*  =  (r*,^*)  is  the  threshold  policy  in  F  with  thresholds  [A*,B*)  defined  by 


(1.46) 


A*  = 


0 


1  -  a 


B*  = 


1-/3 


a 


then 


J{u)  =  Ji(r)  +  J2{6)  >  Ji{r*)  +  =  J(u*). 


(1.47) 


16 


Proof:  This  lemma  follows  from  Corollary  1.2.3  and  Theorem  1.3.1  since 


J2{6)  =  E{C{H,6))  ==  E{C{H,8*))  = 


For  this  problem,  the  requirement  oi{u)  +/3(u)  <  1  is  not  restrictive.  In  fact,  for 
any  uinU  with  Q!(tt)  +  >  1,  the  threshold  policy  u*  =  {r*,6*]  in  F,  defined  by 


r*  0  and 


f  0,  Ci(p  >  C2(l  -  <p) 
\  1,  cicp  <  C2I1  -  ip) 


incurs  lower  cost  than  u. 

The  main  consequence  of  Lemma  1.4.2  is  that  now  in  problem  {Pb),  we  only  have 
to  infimum  over  threshold  policies.  We  now  show  the  infimum  in  (Pb)  is  obtained  by 
a  threshold  policy. 


Theorem  1.4.1.  There  exists  a  threshold  policy  u*  in  F  that  solves  problem  (Pb)- 
The  optimal  thresholds  0<A*<l<B*<oo  with  A*  7^  B*,  are  given  by  the 
relations 


A*  = 


(1.48) 


where  a*  and  b*  are  the  unique  solutions  of  the  transcendental  equations 


C2  +  Ci  =  k(^'(a*)-^'(b*))  (1.49) 

C2(l  -  b*)  =  cia*  +  (b*  -  a*)(ci  -  k^'(a*))  +  k(^(b*)  -  '^(a*)),  (1.50) 

with 

4'(x)  =  (1  —  2a:)  log  — ^ — .  (1-51) 

1  —  a: 

satisfying  0  <  a*  <  b*  <  1. 

Proof;  It  follows  from  Lemma  1.4.1  that  in  order  to  solve  (Pb)  we  only  need  to 
consider  threshold  policies.  Without  loss  in  generality  we  assume  r*  >  0.  From 
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(1.24),  (1.25)  and  (1.42)  the  Bayesian  cost  for  a  threshold  policy  u  is  given  by 
J[u)  =  Ji{t)  +  J2{S) 

ll^flpdi)  +  Ci(p/3[u)  +  C2(l  —  (p)a[u) 
~  ?«(/?(«),  a  (u))  -\-2[\  —  ip)  kw{oL[u),P[u))  +  Cyip  ^[u)  +  C2(l  -  ip)  Oi{u) 

Therefore,  the  optiw.al  threshold  policy  u*  is  given  by 

•^(“*)  •^(«)  (1.52) 

""  kw{^,a)  +  2{l-  <p)  kw{a,0)  +  ci<pl3  +  C2{l-  <p)  a) . 

After  some  algebraic  simplification  using  equations  (1.16),  (1.48),  and  (l.5l), 
equation  (1.52)  becomes 

with 

K{a,  b)  =  {2k^{ip)  +  [(^  -  o(c2(l  -  6)  -  2k^[h)) 

+  {b  -  ip)[cia  -  2k'if{a))  /{b-a)^. 

The  infimum  in  (1.53)  is  over  an  open  set  therefore,  if  the  minimum  value  of  K[a,b) 
exists  then  it  can  be  found  by  solving 

-C2(1-6*)+ci6*)  =0  (1.54) 

-  a*)®' (S’) 

-  C2(l  -  a*)  +  cia*)  =  0  (1.55) 

This  implies  that 


’pEiikJ^  \\ht\\^dt)  +  {l-<p)Eo{kJ^ 


C2  +  Cl  2=  A;('J''(a*)  -  ^''(6*)), 


(1.56) 
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C2(l  -  b*)  =  cia*  +  {b*  -  a*)(ci  -  W{a*))  +  k{<i>{b*)  -  ^(a*)).  (1.57) 

It  can  be  shown  that  for  0<a*<6*<l,  these  equations  have  a  unique  solution 
(Shiryayev  [1977,  pp.  183-184]).  Thus  the  thresholds  are  uniquely  determined  by  ci, 
C2,  and  k.  □ 

1.5  General  Hypopaper  Testing  Problem 

We  will  now  state  the  results  for  the  hypothesis  testing  problem  when  Hq  contains 
a  signal  in  addition  to  noise.  Here,  {xj,  t  >  0}  and  {x|,  t  >  0}  represent  signal 
processes  and  are  and  -valued,  respectively.  The  standard  Brownian  mo¬ 
tions  {wl,  t  >  0}  and  {w|,  t  >  0}  are  R”^^  and  R’^^-valued,  respectively  and  are 
independent  of  the  R^ -valued  standard  Brownian  motion  {u*,  t  >  0}.  Under  each 
hypothesis  the  observed  data,  dyt,  is  the  output  of  a  stochastic  differential  equation 

Under  Hi :  dyt  =  (xj )  dt  +  dvt 

dx\  =  (xj)  dt  -t-  (xj)  dw] 

Under  Hq-.  dyt  =  h^(x|)  dt  -\-  dvt 

dxl^f\xl)dt  +  g\xl)dwl 

The  functions  /^,  /^,  g^,  h^,  and  satisfy  the  Lipschitz  and  growth  conditions 

of  Section  1.1. 

Let  h\  =  Ei{h‘^{x)  \  7/  ),  i  =  1,2.  In  addition  to  (Tl)-(T3)  we  assume  that 

poo 

Pi{  \\hl  -  hlW^  ds  =  oo)  =  1,  ^  =  l,l. 

Jo 

Then,  the  solutions  to  the  Bayesian  and  fixed  probability  of  error  formulations  are 
still  valid  with  A*  defined  by 

A.  =  «xp(  j\hl  -  hi)  dy.  -  \  -  llAjin  ds). 
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1.6  Summary 

We  have  shown  that  in  both  formulations  of  the  sequential  hypothesis  testing  problem, 
the  optimum  decision  policy  is  of  threshold  type.  For  the  fixed  probability  of  error 
formulation,  this  result  was  shown  in  (Liptser  k.  Shiryayev  [1978]).  To  our  knowledge, 
no  one  has  given  explicit  threshold  formulas  for  the  Bayesian  case.  Our  ability  to  do 
so  lies  in  our  modification  of  the  Bayesian  cost.  Usually,  the  Bayesian  cost  includes  a 
constant  penality  for  each  observation.  We  have  chosen  the  observation  cost  to  depend 
on  the  observation  itself.  This  cost  function  is  quite  reasonable  since  it  increases 
the  penality  when  the  confidence  increases  and  results  in  explicit  formulas  for  the 
thresholds. 

In  order  to  implement  the  optimum  policy.  A*  needs  to  be  computed  from  the 
observations  {y^,  s  <  t}.  In  the  next  chapter  we  show  that  under  the  appropriate 
assumptions  on  the  functions  /,  g,  and  h,  the  computation  of  A*  can  be  accomplished 
by  calculating  the  unnormalized  conditional  density  of  x  given  the  observations  and 
then  integrating. 


2.  Numerical  Treatment 


2.1  Introduction 

In  this  chapter  we  will  discuss  the  numerical  method  used  to  approximate  A^.  It  will 
be  shown  that  At  =  u(x,t]  dx  where  u[x,t)  is  the  solution  to  the  Zakai  equation. 
Our  strategy  will  be  to  find  a  good  approximation  to  u(x,  t)  and  then  to  can  use  it  to 
approximate  At.  Therefore  most  of  this  chapter  will  be  spent  approximating  u{x,t). 
Since  the  Zakai  equation  is  a  linear  stochastic  partial  differential  equation,  we  will 
use  standard  numerical  techniques  for  linear  partial  differential  equations  to  obtain 
approximations  to  its  solution.  This  means  using  semigroup  methods  to  prove  the 
necessary  results  and  hence  present  one  convergence  theorem  which  can  be  used  to 
check  convergence  of  several  approximation  schemes. 

This  chapter  is  organized  as  follows:  Section  2  contains  the  necessary  results 
from  nonlinear  filtering.  Sections  3-4  contain  the  necessary  results  from  semigroups 
and  parabolic  partial  differential  equations.  The  remaining  sections  contain  the  ap¬ 
proximation  theorems  for  At  and  u(x,t). 

2.2  Results  from  Nonlinear  Filtering 

From  the  theory  of  nonlinear  filtering,  it  is  known  (Liptser  &;  Shiryayev  [1977])  that 
under  the  appropriate  conditions  on  /,  g,  and  h,  the  unnormalized  density  of  x  given 
the  observations  satisfies  the  linear  stochastic  partial  differential  equation,  known  as 
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the  Zakai  equation,  given  below 

'  du[x,t)  =  L*u[x,t)  dt  +  u[x,t)  h'^ (x)  dyt 
«(x,0)  =  po(a:) 

L*u{x,t)  =  ^  ■^-^[aij{x)u{x,t)\- 

,  (a:,0en 


(2.1) 


where  po{x)  is  the  initial  density  of  x  and  0  =  x  [0,r].  In  order  to  guarantee 
existence  and  uniqueness  of  the  solution  of  (2.1)  (Baras,  Blankenship  &;  Hopkins,  Jr 
[1983]),  we  make  the  following  assumptions  which  are  enforced  throughout. 

(Al)  L*  is  uniformly  elliptic,  that  is,  for  some  A  >  0,  for  all  x,  z  in  R” 


^  cr(x)  z  >  X  z^z. 


(A2)  The  functions  f(x),  cr{x)  and  h[x),  along  with  ^/^(a;), 

■^o'ij{x),  ^hk{x),  a-iid  a^Ta^o’ij(®)  for  hJ  =  l,...,n,  and 

k  =  1, . . .  p,  are  uniformly  bounded  and  Lipschitz  continuous. 

In  order  to  prove  the  necessary  convergence  results  we  will  transform  (2.1)  into 
(2.7)  (given  below),  approximate  the  solution  of  (2.7)  and  then  transform  back  to  get 
an  approximation  to  the  solution  of  (2.1).  This  transformation  is  discussed  in  (Clark 
[1978])  and  is  accomplished  by  a  gauge  transformation  as  follows. 

Under  assumption  (A2)  we  see  that 


L*U  =  ^  (Tij{x) 
»>j=i 


52 


dxidxj 


j=i  ^  i=i  ®  /  3 


where 


j=i  ^3  i=i  ■> 


(2.3) 


To  simplify  the  notation  let 


L*u  =  A*v.  +  c{x)u. 


(2.4) 
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and  define 

<p{3:,yt,t)  =  h'^{x)yt  -  i||/i(x)||^f  +  c(x)t  (2.5) 

Consider  the  function  r[x,t)  which  is  related  to  the  solution  of  (2.1)  via 

r(a:,t)  =  tt(x,t).  (2.6) 


It  follows  from  Ito’s  rule  and  the  definition  of  the  process  {yt,  t  >  0}  that 

r(a:,0)  =  Po(®),  [x^t)  G  fl. 

The  importance  of  this  transformation  is  that  (2.7)  is  a  classical  parabolic  par¬ 
tial  differential  equation  with  coefficients  which  (for  every  oj)  are  Holder  continuous 
functions  of  time.  Therefore  we  can  use  standard  techniques  to  prove  convergence 
of  approximations  of  (2.7)  to  its  solution  and  by  the  simple  transformation  process 
in  (2.6),  get  convergence  of  the  respective  approximation  to  the  solution  of  (2.1).  It 
is  important  to  note  that  the  only  necessary  term  in  the  transformation  is  h'^[x)  yt- 
The  other  terms  in  the  definition  of  (p  help  give  desirable  numerical  properties  to  the 
subsequent  approximation  scheme. 

Once  the  solution  to  (2.1)  is  found,  A*  is  calculated  by  integrating  u[x,t),  i.e., 


(2.8) 


To  see  this  note  that  from  (1.7) 
differential  equation 


and  Ito’s  rule  it  follows  that  At  satisfies  the  stochastic 


dAt  =  AthJ  dyt 


(2.9) 


I,  Ao  —  1. 

Let  <  f,g  >  denote  the  L2-inner  product  of  /  and  g,  i.e., 


<  f,g  >=  /  f{x)g{x)  dx 


(2.10) 


23 


then  from  (2.1) 


d  <  ut,l  >  =<  L*ut,  1  >  dt+  <  Ut  h^,  1  >  dyt 


=<  ut,Ll>  dt+  <  Ut,  h'^  >  dyt 

<ut,h?'  >  , 

=  — - r—  <ut,l>  dyt 

<Ut,l> 

=<  Ut,  1>  hj  dyt. 


(2.11) 


Therefore 


d  <  Ut,l>  —<  Ut,  1  >  h;  dyt 

<  tio,  1  >  =  1- 


(2.12) 


From  existence  and  uniqueness  of  the  solution  of  the  stochastic  differential  equation 
(2.9)  we  obtain 


A-t  — ^  ^ — 


a.s. 


Note  that  in  the  calculation  above,  we  used  the  fact  that  <  Ut,  1  >7^  0.  This  follows 
from  (A1)-(A2). 

In  practical  applications,  xt  and  dyt  correspond  to  physical  signals  and  as  such, 
they  are  bounded.  Therefore  the  boundness  assumptions  in  (A2)  on  /,  g,  and  h 
are  not  restrictive  since  we  can  always  take  them  bounded.  In  fact  because  of  the 
boundness  of  Xt,  we  need  only  consider  values  of  u[x,  t)  for  x  in  a  bounded  set  D  C  M”. 
Therefore,  we  will  solve  (2.1)  with  fi  =  D  x  [0,r]  and  with  the  boundary  condition 
u[x,t)  —  0  for  X  on  dD.  The  assumptions  (Al)  and  (A2)  now  hold  with  replaced 
by  D. 

The  boundness  assumptions  in  (A2)  on  the  partial  derivatives  oi  f ,  cr  ,  and  h  are 
strong.  We  make  these  assumptions  in  order  to  prove  that  each  suitable  approxima¬ 
tion,  Un{x,t),  converges  uniformly  to  the  true  solution  u{x,t)  of  (2.1),  i.e.,  for  each  t 
in  [0,  r]  and  under  the  sup-norm 


lim  |lu„(x,t)  - 'u(x,t)||oo  =  0. 


(2.13) 
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Since  we  will  be  implementing  the  resulting  approximation  scheme  on  a  digital 
computer,  we  must  be  sure  that  under  the  finite  arithmetic  of  the  computer,  A" 
provides  a  good  approximation  to  A*  (with  the  obvious  notation).  Therefore,  we 
insist  on  uniform  convergence.  Notice  that  if  (2.13)  holds  then  by  Holder’s  inequality 
A^  converges  to  At.  We  note  that  under  milder  conditions  on  the  functions  /,  g,  and 
h,  the  corresponding  weak  convergence  results  can  be  shown  (Kushner  [1977]). 
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2.3  Definitions  and  Semigroup  Results 

Let  X  be  a  Banach  space  with  norm  H-H  and  let  Xn  be  a  Banach  space  with  norm 

II  lln' 

Definition  2.3.1.  The  sequence  of  Banach  spaces  approximate  X  if  there 

exist  bounded  linear  operators  Pn  :  X  ^  Xn,  such  that  for  all  f  in  X, 

=  ll/ll-  (2-15) 

Note  that  it  is  not  necessary  for  Xn  to  be  a  subset  of  X.  For  our  purposes  Xn 
will  usually  be  finite  dimensional  while  X  will  always  be  infinite  dimensional. 

Definition  2.3.2.  The  sequence  {/n}o°>  fn  in  Xn,  converges  to  f  in  X,  denoted 

limn,— *00  fn  ~  />  If 

lim  ||/n-Pn/L  =  0.  (2.16) 

Let  B(X,V)  denote  the  set  of  all  bounded  linear  operators  from  X  into  Y  and 
let  B{X)  B{X,X).  The  operator  norm  on  B{X)  will  also  be  denoted  by  ll'H  and 
the  operator  norm  on  B{Xn)  will  be  denoted  H'lljj- 

For  each  t  >  0,  let  T(t)  be  a  continuous  linear  map  in  X. 

Definition  2.3.3.  T(t)  is  said  to  be  a  strongly  continuous  semigroup,  denoted  (Cq)- 
semigroup,  if  it  satisfies  the  following  three  properties: 

(A)  r(o)  =  I 

(B)  T{s  +  t)  =  T{s)-T{t],  s,t>0 

(C)  limtioT{t)f^f  fGX. 

It  is  known  (Yosida  [1980])  that  for  every  (C'o)-semigroup  T(t)  there  exists  con¬ 
stants  M  >  1  and  i  G  ]R  such  that 

||r(t)||  < 


t  >  0. 


(2.17) 
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If  6  =  0  and  M  =  1  we  say  T{t)  is  a  contraction  (Co) -semigroup. 

Definition  2.3.4  The  generator  A  of  a  (Cq) -semigroup  is  defined  to  be  the  unique 
linear  operator  A  satisfying 


lim 

tio 


T{t)f-f 

t 


-Af 


=  0, 


for  all  f  in  D  (A) 


(2.18) 


where  D{A),  the  domain  of  A,  is  defined  by  all  f  in  X  where  the  limit  exists. 

Every  (Co)-semigroup  has  an  infinitesimal  generator,  therefore  we  will  use 
exp(tA)  to  denote  a  (Co)-semigroup  with  infinitesimal  generator  A. 


Lemma  2.3.5.  Let  X  be  a  Banach  space.  Suppose  S  and  S~^  are  bounded  operators 
in  X.  Let  T{t)  be  a  (CQ)-semigroup  in  X  with  infinitesimal  generator  A.  Then 
S~^T{t)  S  is  a  (Co)-semigroup  in  X  with  infinitesimal  generator  S~^AS. 

Proof:  First  we  show  S~^T(t)S  is  a  (Co)-semigroup.  The  conditions  (A)  and  (B) 
being  clearly  satisfied,  we  only  need  to  show  (C)  is  satisfied. 

Let  g  =  Sf  then 

lls-‘r(i)  Sf-f\\<  lls-‘  II  ||T(()  Sf  -  s/ll 

=  l|s-‘ll  llrWs-sHo- 


Therefore,  S  ^T{t)  S  is  a  (Co)-semigroup,  and  we  now  show  that  S  ^AS  is  its  in¬ 
finitesimal  generator.  For  all  /  in  X  such  that  S/eD  (A) 


lim 


S-'^T{t)Sf-  f 


-  {S-'^AS)f 


lim 

tio 


^_,fT{t)Sf-Sf 


AS 


f) 


=  11-5~^||  lim 


T{t)g-g 


-Ag 


0. 


□ 


Definition  2.3.6.  A  linear  operator  A  is  closed  if  the  graph  of  A  dehned  by 


^(A)  =  {(/,A/):  fED(A)}  (2.19) 


is  a  closed  subspace  of  X  X  X. 
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Definition  2.3.7.  A  linear  operator  A  in  a  Banach  space  X  is  closeable  if  the  closure 
of  the  graph  Q  (A)  is  itself  the  graph  of  a  linear  operator  in  X.  The  closure  of  A  is 
denoted  A. 

Definition  2.3.8.  Let  the  sequence  of  Banach  spaces  approximate  X.  The 

linear  operators  A„  in  Xn  converge  to  the  linear  operator  A  in  X,  denoted  A  = 
lim„-H.oo  An,  if  for  all  f  in  D(A) 

lim  \\AnPnf-PnAf\\^  =  0.  (2.20) 


The  following  theorem  from  Hille  and  Yosida  characterizes  the  linear  operators 
in  X  which  generate  (Co)-semigroups. 

Theorem  2.3.9  (Hille- Yosida).  Let  A  be  a  closed  linear  operator  with  domain  and 
range  in  a  Banach  space  X.  Let  b  E  R  and  M  >  1.  The  following  assertions  are 
equivalent: 

(i)  The  operator  A  generates  a  (Co) -semigroup,  T(t)  for  which 

||r(t)||  <  Me'’*  t>0.  (2.21) 

(ii)  The  domain  of  A  is  dense  in  X,  sI-A  is  boundedly  invertible  for  s  >  b  and 

||(s-6)”(s/- A)“"||  <  M,  nelX  (2.22) 

Proof:  See  (Yosida  [1980,  pp.  246  -248  ]). 

The  following  approximation  theorem  for  (Cq) -semigroups  comes  from  Trotter. 

Theorem  2.3.10  (Trotter).  Let  {hn}'o  be  a  sequence  of  positive  numbers  converging 
to  zero,  and  {T'n}o°  ®  sequence  of  linear  operators  in  Xn  satisfying  the  stability 


condition 


(2.23) 
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where  M  and  b  are  constants,  independent  of  n  and  k.  Let  An  =  h~^{Tn  —  I)  so  that 
r*  =  (/  +  hnAn)^  and  define  A  :=  lim„_Koo  An-  If 
(i)  P{A)  is  dense  in  X,  and 

(a)  For  some  s  >  b,  the  range  of  (si  —  A)  is  dense  in  X, 
then  the  closure  of  A  is  the  inhnitesimal  generator  of  a  (Co)-semigroup  T(t)  and 

r(i)  =  lim  t>0.  (2.24) 

n—^oo 


Proof:  See  (Trotter  [1958,  pp.  903-904  ]). 


Since  the  operator  A(t)  in  (2.7)  depends  on  time  we  are  interested  in  approximat¬ 
ing  the  solution  of  the  abstract  Cauchy  problem  for  linear  problems.  In  the  Cauchy 
problem  we  want  to  find  u(t)  in  X  such  that 


■u(O)  =  uo 


(2.26) 


where  for  each  t  >  0,  A(t)  is  a  linear  operator  in  a  Banach  space  X.  We  assume  that 
(2.26)  is  well-posed,  i.e.,  there  exists  a  unique  solution  to  (2.26)  in  X. 

If  A{t)  =  A  then  the  solution  u(t)  to  (2.26)  is  given  by 


u(t)  =  T{t)  Uo 


(2.27) 


where  A  is  the  infinitesimal  generator  of  the  semigroup  T{t).  In  this  case.  Theo¬ 
rem  2.3.10  describes  how  to  approximate  T(t)  and  hence  how  to  approximate  u(t). 

Definition  2.3.12.  The  operator  U(t,s)  is  called  a  fundamental  solution  to  (2.26)  if 
it  satisGes 

(1)  U(t,s)  is  a  strongly  continuous  function,  deGned  for  0  <  s  <  t  <  T  which 
takes  values  in  B(X). 

(2)  U  (t,  r)  U  (r,  s)  =  U  (t,  s)  for  0<s<r<t<T. 
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(3)  U(s,s)  =  /  foraUse  lO.T). 

W  T,U{t,s)=A(t)U{t,s). 

(5)  £U{t,B)  =  -U{t,s)A(s). 

Conditions  (4)-(5)  are  understood  to  hold  on  a  dense  subspace  of  X  where  they 
make  sense.  The  derivatives  ^  and  ^  are  taken  in  the  strong  topology  of  X.  Notice 
that  if  A{t)  =  A  then  U[t,s)  is  just  the  (Co) -semigroup  T[t  —  s)  with  generator  A. 
The  solution  to  (2.26)  is  given  by 

u{t)  =U{t,0)uo.  (2.28) 

where  U{t,s)  is  the  fundamental  solution  to  (2.26). 

Let  X  be  a  Banach  space  with  norm  H-H.  Let  T  be  a  dense  subspace  of  X  and 
assume  Y  is  itself  a  Banach  space,  with  norm  ||•||y.  Suppose  there  exists  a  constant  c 
such  that  ||u||  <  c||u||j,  for  all  v  in  Y.  Henceforth,  we  assume  X  and  Y  satisfy  these 
conditions.  We  now  give  some  useful  definitions. 

Definition  2.3.13.  The  set  of  all  generators  of  (Co)-semigroups  with  constants  M 
and  b  in  X  is  denoted  by  by  G{X,M,b)  ,  and  the  set  of  all  generators  of  (Cq)- 
semigroups  in  X  is  denoted  by 

U  U  G(X,M,b)  (2.29) 

6gIRM>1 


Definition  2.3.14.  Let  A  G  C(X).  The  Banach  space  Y  is  called  A-admissible  if 
exp{tA)  maps  Y  into  Y,  and  if  the  restriction  of  exp(t  A)  to  Y  forms  a  semigroup 
in  Y. 

Definition  2.3.15.  The  family  of  operators  A{t)  in  G{X)  is  called  stable,  with 
stability  constants  M  and  b,  if  there  exist  real  numbers  M  >  1  and  b  such  that 


k 

n  exp(sj  A(tj)) 
j=i 


Sj  >  0, 


(2.30) 
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for  all  0  <  t\  <  t2  <  •  •  •  <  tk  <T  and  k  —  1,2,....  Here  the  product  is  time  ordered, 

i.e. 

k 

JJ  exp{sj  A{tj) )  =  exp{sk  A{tk) )  •  •  •  exp( si  A{ti) ).  (2.31) 

3-i 

For  the  remainder  of  this  chapter,  all  operator  products  will  be  time  ordered. 
The  following  lemma  will  be  useful  to  test  for  stability. 

Lemma  2.3.16.  Assume  that  A[t)  is  stable  with  stability  constants  M  and  b.  If 
for  each  t  in  [0,  T],  B(t)  is  bounded,  i.e.,  ||B(t)||  <  K  <  oo  then  for  each  t  in  [0,  i], 
A(t)  +  B(t)  belongs  to  G(X)  and  is  stable  with  stability  constants  M  and  b  +  MK, 

Proof:  See  (Kato  [1970,  p  248  ]). 

We  now  state  assumptions  on  the  generator  A{t)  in  (2.26). 

(A3)  A{t)  is  stable  with  constants  1  and  b. 

(A4)  The  Banach  space  Y  is  A-admissible  and  Y  C  V[A[t))  for  each  t  in  [0,  T] 
and  the  operator  A[t)  is  a  continuous  function  in  the  norm  of  B(Y,X). 

(A5)  A(t)  is  closed  in  X. 

(A6)  For  some  s  >  b,  the  range  of  {sI—A{t))  is  dense  in  X,  where  b  is  independent 
of  t. 

The  following  convergence  theorem  combines  the  theorems  (Kato  [l970,Theorem 
4.1  p247])  and  (Trotter  [1958,  Theorem  5.3  p903]). 

Theorem  2.3.17.  Let  hn  be  a  sequence  of  positive  numbers  such  that  hn  0  as 
n  — >  oo.  Suppose  that  for  each  t  in  [O,  T]  the  linear  operators  Tn,t  ia  satisfy  the 
stability  conditions 

(2.32a) 

and 

N 

WTr^.tj  (2.32b) 

3=0 
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where  tj  =  j  hn,  <T  and  the  constants  M  >  1,  b  E  IR  are  independent  of  t,  k  and 

n. 

Let  An{t)  =  {Tn,t  —  I)/ hn  and  suppose  that 


A{t)  :=  lim  An{t). 


(2.33) 


for  each  tk  in  [0,r].  For  each  pair  {t,s)  satisfying  0<s<t<T,tk<s<  tk+i,  and 
ti  <  t  <  define  the  operator 


(I  A:  =  /, 


C/n(i,a)  =  < 


i-i 

n  <  '• 


(2.34) 


I  3=k 

If  the  operator  A{t)  in  (2.33)  satisfies  (A3)-(A6)  and  if  the  corresponding  Cauchy 
problem  in  well-posed  then 


U(t,s)=  lim  Un(t,s)  (2.35) 

where  U(t,  s)  is  the  fundamental  solution  generated  by  A{t). 

Proof:  It  is  easy  to  see  from  its  definition  that  Un[t,s)  satisfies  (l),  (2),  and  (3)  of 
Definition  2.3.12.  Furthermore  it  can  be  seen  that  t7„(t,  s)  maps  Xn  into  Xn,  and 
that  for  each  n  in  F, 

M  +  '‘2^)  -  =  A4t)  y„(f,  s)P„v  (2.37) 

and 

hn) - U{t,  ^  -Un{t,  S  +  hn)  An{s)PnV.  (2.38) 

hn 

From  (2.32)  and  the  stability  of  A{t)  we  have  that 

\\Un{t,s)\\<Me^^*-^\  ||C/(f,s)l|^  (2.39) 

Next  V  eY  and  let 

U{t,s)  =U{hn[t/hn],hn[s/hn\). 
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It  follows  that  all  we  need  to  show  is  the  convergence  of  {t,  s)  to  U  (t,  s)  hence 


l-i 

\\Un{t,s)PnV  -  Fni7(f,s)t;||^  =  \\Y2Un{t,tj+i)PnU{tj  +  i,s)v  -  Un{t,tj)PnU{tj,s)v\\^ 

j=k 

1-1 

+  “  Un{t,tj  +  l)PnU{tj,s)v 

j=k 

+  Ur,{t,tj+i)PnU{tj,s)v  -  Un{t,tj)PnU{tj,s)v\\^ 

1-1 

=  \\'^Ur,{t,tj+i)Pn{U{tj  +  l,s)  -U{tj,s))v 

}=k 

+  (Unit.tj+i)  -  Ur,{t,tj))PnUitj,s)v\\^ 

1-1 

=  \\^Un{t,tj^l)Pn{U{tj  +  utj)  -  I)U[tj,s)v 

j=k 

—  Unit,  tj^i)hnAn{tj)PnUitj,  s)t;  ||^ 

1-1 

<  -  I>  -  hnAn{tj)Pnv\\^ 

j=k 

3„p  (2.40) 

j=k,...,l-l  hn 

where  7  =  max{6,6}  and  M'  =  MM.  From  (2.33)  and  Definition  2.3.12  it  follows 
that  Unit,  s)PnV  converges  as  n.  00  uniformly  in  0  <  s  <  t  <  T".  From  the  fact  that 
Y  is  dense  in  X  it  follows  that  for  all  u  in  X, 


lim  \\Unit,s)PnU-PnU{t,s)u\\  =  0  (2.4l) 

n—^oo 

exists  uniformly  in0<s<t<r.  □ 

The  results  above  can  be  used  to  show  convergence  of  several  different  types  of 
approximation  schemes.  In  the  next  sections,  we  will  apply  them  to  implicit  Euler 
type,  finite  difference  approximations  for  the  parabolic  equation  (2.7). 
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2.4  Results  on  Parabolic  Partial  Differential  Equations 

In  this  section  we  will  quote  results  from  the  theory  of  parabolic  partial  differential 
equations  to  establish  that  A*  in  (2.4)  generates  a  (C'o)-semigroup  and  to  show  that 
there  exists  a  unique  solution  to  (2.7).  This  will  done  for  the  case  where  x  in  (2.7)  is 
a  scalar  diffusion.  Similar  techniques  can  be  applied  to  the  vector  case  and  results  in 
schemes  similar  to  those  found  in  (Richtmyer  &;  Morton  [1967]). 

Throughout  this  section,  let  D  =  (a,  b)  a  bounded  interval  and  let  X  =  Ch[[a,  6)]n 
Co  [(a,  6)1  under  the  sup-norm  and  Y  =  C°°[(a,  6)]  Pi  Co  [(a,  6)]  with  norm  |[■||y  such 
that  for  each  /  in  Y 

OO 

ll/lb  =  E  IIP' II- 

n=0 

It  is  known  that  Y  forms  a  dense  subspace  in  X.  Furthermore,  for  every  /  in  Y, 
1|/||  <  l|/||y  We  will  consider  the  Cauchy  problem: 

'  ut{x,t)  =  a{x,t)  Uxx{x,t)  +  b{x,t)  Ux{x,t)  +  c{x,t)  u{x,t), 

<  u{a,t)  =  u{b,t)  =  0,  t€[0,r]  (2.43) 

,  ti(x,  O)  =  uo(a:),  (x,t)  G  fl, 

where  fl  =  (a,  b)  X  [0,  T].  The  follow  theorem  from  Besala  gives  the  necessary  existence 
result  for  the  solution  of  (2.43). 

Theorem  2.4.1  (Besala)  For  (2.43),  let  the  functions  a{x,t),  b{x,t),  and  c(x,t)  (real 
valued)  together  with  ax{x,t),  aa;x(x,t),  bx{x,t)  be  locally  Holder  continuous  in  U. 
Assume  that  for  all  {x,  t)  in  fl 

i)  a{x,t)  >  A  >  0,  for  some  X 
a)  c(x,t)  <  0, 

Hi)  c(x,t)  -  bx{x,t)  +axx{x,t)  <  0. 

Then  the  Cauchy  problem  (2.43)  has  a  fundamental  solution  U  (t,  s)  in  the  Banach 


space  X  where 


||t/(t,s)||<l. 


(2.44) 
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Furthermore,  there  exists  a  function  z,  s)  which  satisfies 


0  <  r(a:,  s)  < 


-  s) 


(2.45) 


for  some  positive  k,  and 


/  T{x,t;z,s)  dz  <  1 
J  a 

fb 

/  s)  dx  <  1 

J  a 


such  that 


U{t,s)f{x)=j  T{x,t-,z,s)  f{z)  dz. 

J  a 

Moreover,  if  uo(x)  is  continuous  and  bounded,  then 

u[x,t)  =  U{t,0)uo{x)  =  I  T  {x,t\  z,  0)  uq  {z)  dz 

J  a 


(2.46) 


(2.47) 


is  a  bounded  solution  of  (2.43). 

Proof:  See  (Besala  [1975]). 

The  following  theorem  due  to  Friedman  gives  the  necessary  uniqueness  result  for 
(2.43). 

Theorem  2.4.2.  Ifa{x,t),  b{x,t),  and  c(x,t),  in  (2.43)  are  continuous  in  ft.  If  there 
exists  A  >  0  such  that  for  all  x  in  (a,  b) 


a{x,t)  >  X  (2.48) 

for  {x,t)  in  fi,  then  there  exists  at  most  one  solution  to  the  Cauchy  problem  (2.43). 
Proof:  See  (Friedman  [1964,  Theorem  7,  p.  41]) 

Corollary  2.4.3.  For  each  tk  in  [0,r],  the  operator  A[tk)  defined  in  (2.7)  generates 
a  (Co )-semigroup  in  X  with  constants  1  and  K. 
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Proof:  Using  the  notation  of  (2.43)  for  the  operator  A[tk)  and  from  assumption  (A2) 
it  follows  that  there  exists  a  finite  constant  K  >  0  such  that 


m.a.x{axx(x,t)  —  bx{x,t)  —  K}  <  0.  (2.50) 

aelR 

Define  A[tk,)  =  A[tk)  —  K.  From  assumption  (Al)  it  follows  that  A[tk)  satisfies  the 
hypotheses  of  Theorem  2.4.1,  therefore  A{tk)  generates  a  contraction  (Co)-semigroup. 
From  Lemma  2.3.16  it  follows  that  that  A{tk)  generates  a  (C'o)-semigroup  in  X  with 
constants  1  and  K. 

Corollary  2.4.4.  Under  assumptions  (Al)-(A2)  the  semigroup  generated  by  A(tk) 
for  each  tk  in  [0,r]  is  given  by 

exp(t  A(tk) )  =  exp(t  A* )  (2-51) 


where  is  defined  in  (2.23). 

Proof:  Directly  from  Lemma  2.3.5. 

We  are  now  ready  to  show  a  convergence  theorem  for  the  case  when  only  time  is 
discretized 


Theorem  2.4.5.  Let  At  =  T/n  and  k  =  [t/At].  For  k  >  0  define 


Un{t,0)poix)  =  JJ  \ 

^  j=o  ' 

(2.52) 

where  Ayy  =  y{j+i)At  -  VjAt,  then 


k-i 


lim  sup 

n->cx3  tg[o,T] 


Un{t,0)po{x)  -  r{x,t) 


(2.53) 


where  r{x,t)  is  the  solution  of  (2.7).  Furthermore,  in  view  of  (2.6),  it  follows  that 


lim  sup 
n-*oo  te[o,r] 


e'^(='’^'=^*’^^‘)£7„(i,0)po(a:)  -  u{x,t) 
where  u(x,t)  is  the  solution  of  (2.1). 


(2.54) 
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Proof:  Let  Tn,t  =  exp(AtA(t))  where  A{t)  is  the  operator  in  (2.7).  Substituting 
(2.51)  into  the  definition  of  I7„(t, 0)  of  Corollary  2.3.18  and  rewriting  the  product 
yields  (2.52).  From  Corollary  2.3.4,  A[t)  is  stable  with  constants  1  and  K  hence  all 
of  the  hypothesis  of  Theorem  2.3.17  are  satisfied,  therefore  the  convergence  (2.53). 


Corollary  2.4.6.  Let  An{t)  A(tt+i)  for  tk  <  t  <  tk+i  replace  An(t)  in  Theo¬ 
rem  2.3.17.  Iftk  =  k  At  then  for  A(t)  in  (2.7)  and  using  (2.6),  we  see  that 


V4t,0)poix)  =  (  n  po(a:)  (2.55) 

i=o 


converges  to  u[x,t)  in  (2.1),  i.e.. 


lim  sup  ||Pn(i,0)po(a:)  -ti(x,t)||  =  0. 

t6[0,T] 


(2.56) 


Proof:  All  of  the  hypotheses  of  Theorem  2.3.17  are  satisfied  with  An{t),  hence  the 
result  (2.35)  with  only  the  indices  changed. 


Example  1.  Time  discretization  of  (2.1). 

Under  the  assumptions  (Al)-(A2)  the  discrete  time  approximation  scheme 

Vn{x,t)=  (  (2.57) 

^3=0  ' 


converges  to  the  solution  of  (2.1). 

Let  Tn  =  {I—AtL*)^^  and  the  Banach  space  X  as  above.  It  is  known  (Richtmyer 
h  Morton  [1967])  that 

llTnli  <  (2.58) 


therefore 

exp(tL*)  =  lim  (r„)[^]. 


(2.59) 
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Convergence  of  (2.57)  to  the  solution  of  (2.1)  follows  from  Theorem  2.3.17  and  Corol¬ 
lary  2.5.6.  This  scheme  appears  in  (LeGland  [1981];  Pardoux  &  Talay  [1983])  and 
was  obtained  by  different  methods. 

By  arguments  similar  to  those  in  Theorem  2.4.5  it  is  possible  to  show  that 


Ki(t,0)p(a:)  = 


)po(a:) 


(2.60) 


converges  to  the  solution  of  (2.1),  where  Aj/y  and  At  are  defined  in  Theorem  2.4.5. 
This  result  is  important  for  applications  since  it  demonstrates  the  importance  of  the 
transition  density  of  x  in  the  calculation  of  u{x,t).  The  importance  of  the  transition 
density  has  been  shown  and  is  illustrated  in  the  Kallianpur-Striebel  formula.  For  us,  it 
is  interesting  that  (2.60)  can  be  thought  of  as  a  direct  approximation  the  Kallianpur- 
Striebel  formula. 


In  applications,  the  transition  density  of  x  can  usually  be  approximated  from 
the  observed  data.  While  the  functions  /  and  g  usually  cannot.  Therefore,  this 
approximation  method  has  the  additional,  practical  advantage  of  just  needing  the 
transition  density. 
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2.5  The  Finite  Difference  Approximation  Scheme  for  (2.1) 


It  is  now  possible  to  find  several  finite  difference  approximation  schemes  for  (2.1). 
In  fact,  from  Theorem  2.4.5  it  follows  that  all  we  need  is  a  good  approximation  for 
exp(AtA*).  Several  good  approximation  schemes  exist  in  the  numerical  analysis 
literature,  however  we  will  only  concentrate  on  one. 

In  the  scalar  x  case  we  have  D  =  (a,  b) , 

L*u[x,t)  =  a{x)  Uxx{x,t)  +  b[x)  Ux{x,t)  +  c[x)u[x,t) 

(2.60) 


=  A* u{x,  t)  +  c{x)u{x,  t) 


where 

b{x)  =  g{x)  g'{x)  -  f{x)  (2.61) 


c{x)  =  ^(x)  g"{x)  +  (ff'(x))^  +  g{x)  g'{x)  -  f'{x) 

Let  Ax  >  0  and  define  x*,  =  c  +  kAx  and  n  such  that  x^  <  b.  Consider  the 

collection  of  points  {xfc}o  in  D.  Let  Ax  0  as  n  ^  oo  such  that  Xn  ^  b  as  n  —>■  oo. 
Let  X  and  Y  be  the  Banach  spaces  in  Section  2.4  and  let  the  approximating  Banach 
spaces  Xn  =  under  the  oo-norm.  Define  the  operator  Pn  :  X  Xn  such  that 

for  each  in  A 

{Pn<t>)i  =  ^[xi),  i  =  0,...,n.  (2.62) 


It  is  easily  checked  that  for  each  in  X 


lim  ||Pn<?^||^  =  . 


The  linear  operator  A„  in  Xn  will  be  obtained  from  A*  by  replacing  the  x- 
derivatives  in  A*  with  finite  difference  approximations.  To  this  end,  for  each  (f>  in 
Y  pose 

,  ,  J\  r  ^(f){xi+i) -2<f>{xi) +<f>{xi-l)  _  (l){xi+i)  -  <j){xi) 

{AnPn<f>)i  =  a{xi)^ - ^ - h  max(6(xi) ,  0) 


(Ax): 


+  min(6(xj),0) 


Ax 

<f){xi)  -  <f){xi-l) 
Ax 


(2.63) 
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where  i=0,. . .  ,n. 

From  Taylor’s  Theorem 

<f>{xi  ±  Aa:)  =  <t>{xi)  +  <;6a;(a:i)(±Aa:)  +  )-cj)xx{xi  +  ^±)(Aa:)^  (2.64) 

z 

where  0  <  <  Aa:  and  —  Aa:  <  <  0.  Substituting  into  (2.63)  yields 

{AnPn(l))i  =  {(j>xx{xi  +  ^4-)  +  <i>xx{xi  +  OJ))  +  b{xi)  {<j)x{xi)  +  0(Aa:)).  (2.65) 

Since  <f>xx  is  continuous  it  follows  that 


lim  \\Ar,P,,<j>-PnA*<l>\\^  =  0. 


(2.66) 


(2.67) 


We  are  now  able  to  show  convergence  of  the  full  discretization  of  (2.7),  i.e.  when 
both  time  and  space  are  discretized. 

Theorem  2.5.1.  Let  At  =  hn  and  let  Pn  as  in  (2.62).  Define 
Dk  =  diag 

=  diag  ^yk+{e(.xi)-^h^{xi)) 

where  Ayk  =  y{k+i)At  ~  VkAt-  Consider  v*  in  defined  by 

yfc+i  ^  (/_ 

V°  =  PnPO- 

Let  kAt  — >  t  as  n  — >  oo,  then  for  every  t  in  [O,  T] 


(2.68) 


lim  sup  |(v*^)i  —  u(xi,  ^At)  I  —  0. 

n-^oo  ,=o,...,n 


(2.69) 


Proof:  Let  Tn  =  {I  —  At  An)  From  (2.63)  we  see  that  An  is  a  diagonally  dominant 
matrix.  It  has  the  form 

/-  +  A 


A 


n  — 


+ 


(2.70) 


V 


+  -y 
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For  every  Ai  >  0,  {Tn)  ^  has  the  form 


'+  - 


(r„)-‘  =  /-A(A„  = 


(2.71) 


and  is  strictly  diagonally  dominant  since  the  diagonal  is  reinforced.  It  is  easily  checked 
that  is  an  inverse  positive  matrix,  i.e.,  {Tn)ij  >  0,  see  (Schroder  [1978,  Corol¬ 
lary  1.6b,  p221  ]).  Since  is  invertible  for  all  At  in  [0, T]  we  see  that 

||T„|U  < 

for  some  b  in  IR,  so  that  Tn  satisfies  (2.23). 

Let  ~  [Ek]~^  where 


FJfc  =  diag 
=  diag 

and  p(',-,-)  is  defined  in  (2.5).  Then  from  (2.68)  we  obtain 

I  «;°  =  po 

Define  Tn{tk)  —  {Ek+\)~^  Tn  S'fc+i  where  tk  =  kAt  then  for  each  tk  in  [0,T] 


(2.72) 


(2.73) 


l[Tn{tj)  <Me‘ 


As  in  Theorem  2.3.10,  we  define 


An{tk)  = 


Tn{tk)  -  I 


(2.74) 
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From  the  continuity  and  differentiability  of  <p{-,  •,  •)  and  from  (2.66)  it  follows  that  for 
each  <!>  in  D{A*) 

n— »>oo  '  ^  M  oo 

Therefore  from  Theorem  2.3.17,  in  (2.73)  converges  to  the  solution  of  (2.7).  The 
continuity  of  the  transformation  (2.6)  implies  (2.68)  converges  to  the  solution  of  (2.1) 
hence  (2.69).  □ 

The  approximation  for  b{x)  u(x,t)  is  standard  in  the  numerical  literature  and  is 
motivated  from  stability  of  finite  difference  methods  for  hyperbolic  equations  (Richt- 
myer  Ac  Morton  [1967,  p.  292  ]).  Basically,  when  b[xi)  >  0  we  take  a  forward  difference 
approximation  and  when  b[xi)  <  0  we  take  a  backward  difference  approximation.  Us¬ 
ing  this  approximation  ensures  positivity  of  our  approximation  scheme,  independent 
of  the  value  of  At. 
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2.6  Convergence  of  the  Corresponding  Likelihood  Ratio  Approximation 

Using  defined  in  (2.69)  it  is  easy  to  construct  a  convergent  approximation  to  the 
likelihood  ratio.  Let  be  the  collection  of  points  defined  in  Section  2.5  and  let 


Un(x,t)  =  (n*)i  if  Xi  <  X  <  Si+i,  At  <  t  <  (A:  +  l)At 


and  define 


A? 


n 

dx  = 

i=0 


then 


lim  |At  —  A” 


<  (b  —  a)  lim  sup  \u{x,i)  —  Un{x,t)\ 

=  (b  —  a)  lim  sup  \u[xi,t)  —  Un{xi,t)\ 

=  (b  —  a)  lim  sup  \u{xi,t)  -  {v’^)i\ 

i=0 . n 


=  0. 


Therefore,  A”  is  a  convergent  approximation  to  A^. 
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2.7  Summary 

The  main  result  in  chapter  was  the  convergence  result  in  Theorem  2.3.17.  This 
theorem  presented  verifiable  conditions  which  ensure  convergence  of  an  approximation 
scheme  for  the  abstract  Cauchy  problem  (2.26).  These  conditions  were  used  to  show 
convergence  of  the  natural  finite  difference  approximation  scheme  for  the  solution  to 
(2.1).  This,  in  turn,  implied  show  convergence  of  the  approximate  likelihood  ratio  to 
the  true  likelihood  ratio. 

Schemes  similar  to  (2.68)  have  been  discussed  in  (Kushner  [1977];  Pardoux  h 
Talay  [1983]).  Furthermore,  numerical  studies  have  been  performed  in  (Yavin  [1985]) 
using  these  methods  which  have  produced  satisfactory  results  for  approximations  to 

A 

ht. 

Again,  the  importance  of  the  result  in  Theorem  2.4.5  is  that  the  stochastic  part 
of  (2.1)  has  been  isolated  from  the  non-stochastic  part.  Therefore,  the  wealth  of 
information  on  numerical  approximation  for  parabolic  partial  differential  equations 
can  be  used  to  approximate  exp(tA*  ),  and  this  in  turn  is  used  to  approximate  (2.1). 


3.  VLSI  Architectures 


3.1  Introduction 

In  this  chapter  we  will  present  a  VLSI  architecture  for  solving  (2.1)  when  x  is  scalar. 
This  architecture  will  not  be  efficient  for  vector  x.  However,  as  was  pointed  out  in 
Section  2.9,  it  is  possible  to  use  some  of  the  advanced  techniques  for  sparse  matrices  to 
develop  architectures  for  higher  state  dimensions.  It  is  important  to  emphasize  that 
these  methods  are  not  sufficiently  efficient  for  real  time  implementations  of  problems 
with  state  dimensions  higher  than  three. 

In  Section  2.7  we  presented  a  finite  difference  scheme  to  approximate  the  solution 
of  (2.1).  At  each  time  t  —  kAt,  this  scheme  involved  solving  the  linear  equation 

(/- A<A,,)V^+^  =  L>feV^  (3.1) 

where  Dk  =  diag{exp(/i(a:i)  A2/fc  +  (c(a:i)  -  ^h^{xi))  At)}.  Our  goal  is  to  give  a  design 
for  a  VLSI  chip  to  efficiently  solve  (3.1).  This  means  that: 

(1)  the  time  necessary  to  compute  given  V*,  A,  and  yk,  should  be  below 

a  problem  dependent  threshold;  and 

(2)  the  control  structure  within  the  chip  should  be  simple  and  regular. 

We  will  show  how  the  systolic  array  architecture  of  (Kung  &  Leiserson  [1980]) 
can  be  used  to  obtain  a  VLSI  design  satisfying  these  goals.  The  discussion  of  systolic 
architectures  follows  closely  (Kung  Sc  Leiserson  [1980]). 
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3.2  Systolic  Arrays  and  Sequential  Computations 


Systolic  processors  are  arrays  of  tightly  synchronized  simple  processors  in  which  data 
is  fed  to  each  processor  in  a  regular  and  ordered  manner.  In  these  arrays,  data  flows 
through  the  processors  much  in  the  same  way  that  blood  flows  through  the  heart, 
hence  the  term  systolic.  Because  of  their  order  and  regularity,  systolic  processors  are 
significantly  faster  than  conventional  processors. 

We  will  be  interested  in  systolic  processors  which  perform  linear  algebra  oper¬ 
ations.  The  basic  component  of  these  arrays  is  the  inner  product  processor  (IPP) 
X  y 


denoted  by 


T 


.  At  each  clock  pulse  the  IPP  takes  the  inputs  x,  y  and  a  and 


computes  ax-\-y  (the  inner  product  step).  This  value  is  output  on  the  y-output  line, 
and  the  x  and  a  values  pass  through  to  their  respective  output  lines  untouched. 

To  illustrate  the  speedup,  consider  how  a  typical  von  Neumann  computer  would 
perform  the  inner  product  step.  Before  a  von  Neumann  computer  can  execute  an 
instruction  it  must  first  get  (or  fetch)  the  instruction  from  memory  and  then  fetch 
the  arguments  (or  operands)  for  that  instruction  from  memory.  Table  3.2.1  shows 
the  computer  instructions  necessary  to  execute  the  inner  product  step.  It  does  not 
include  the  overhead  of  fetching  the  instructions  and  the  operands. 


Instruction 

(  cont.  ) 

1.  Fetch  X 

1.  Fetch  a 

2.  Multiply  a  X 

3.  Fetch  y 

5.  Add  (ca:)  ii  y 

6.  return  {ax  -f  y) 

Table  3.2.1.  Operations  performed  executing  the  inner  product  step. 


However,  the  computer  actually  performs  two  additional  operations  in  executing 
each  imand.  For  each  instruction  listed  in  Table  3.2.1  .  t:  are  three  additional 
instructions  to  perform.  Table  3.2.2  shows  these  additional  instructions. 

Therefore,  performing  the  inner  product  step  requires  a  von  Neumann  computer 
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Actual  Operations 

Fetch  Instruction 

Fetch  Operand  for  Instruction 

Execute  Instruction 


Table  3.2.2.  Operations  performed  each  instruction  in  Table  S.l. 


to  execute  18  operations.  Notice  that  only  two  of  these  operations  are  mathematical, 
the  other  16  involve  data  manipulation  and  program  control.  The  time  required  to 
perform  each  of  the  18  operations  varies  from  computer  to  computer,  however  on  any 
one  computer,  these  steps  take  approximately  the  same  time  to  perform  (when  using 
fixed  point  arithmatic).  Systolic  processors  are  fast  because  they  eliminate  the  16 
wasted  operations. 

To  illustrate  how  a  systolic  array  can  be  used  in  a  linear  algebra  operation  we 
consider  the  matrix-vector  operation  y  =  Ax  where  A  is  an  n-by-n  matrix  and  a:  is  a 
vector  in  Let  pi  denote  the  i*^  component  of  the  vector  y.  The  following  n-step 
recursion  can  be  used  to  calculate  yi 

=  0 

=  yj*"^  +  aik  xk,  l<k<n. 

It  is  easily  verified  that  y^  = 

The  above  recursion  shows  that  matrix-vector  multiplication  is  nothing  more 
than  a  sum  of  inner  product  steps.  Therefore,  it  is  reasonable  to  expect  that  with  the 
correct  data  flow  we  can  accomplish  this  matrix-vector  multiplication  with  an  array 
of  IPP’s. 

Let  A  be  an  n-by-n  matrix.  We  say  A  is  a  band  matrix  with  bandwidth  w  if 
aij  =  0  outside  a  w-wide  diagonal  band  containing  the  main  diagonal.  For  example. 
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the  matrix  A  in  (3.2)  has  bandwidth  3. 


/ai.i 

“1,2 

0 

0 

“2,1 

“2,2 

“2,3 

0 

0 

“3.1 

“3,2 

“3,3 

“3,4 

0 

0 

“4,2 

“4,3 

“4,4 

“4,6 

V  0 

0 
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0 
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0 

0 

^33 

0 

^42 

“23 

0 

“32 

0 

0 

a22 

0 

<^31 

OI2 

0 

021 

0 

0 

Oil 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

(3.2) 


yi  o  ya  °  ya  •  •  • 
X  output 


Figure  3.2.3.  Systolic  layout  for  matrix  vector  multiplication. 


Figure  3.2.3  from  (Kung  &  Leiserson  [1980])  shows  how  to  compute  y  —  Ax  for 
A  in  (3.2).  Each  j/*  is  initially  zero,  and  stands  for  a  one  step  delay.  Figure  3.2.4 
from  (Kung  &  Leiserson  [1980]),  shows  the  array  at  several  clock  steps  to  show  how 
the  data  flow  and  interconnections  in  Figure  3.2.3  calculate  y. 

For  our  example,  A  had  bandwidth  3.  In  general,  if  A  has  bandwidth  w  then 
we  need  w  IPPs  to  accomplish  a  matrix-vector  multiplication.  It  takes  2n  •+■  tu  clock 
cycles  for  this  algorithm  to  compute  the  matrix  vector  multiplication  as  compared 
to  the  sequential  algorithm  which  takes  0{nw)  units  of  time.  As  pointed  out  by 
Tables  3.2. 1-3. 2.2,  a  unit  of  time  for  the  sequential  algorithm  is  several  clock  cycles 
long.  Therefore,  the  systolic  algorithm  is  significantly  faster. 

The  number  of  processors  depends  only  on  the  bandwidth  of  the  matrix  A  and 
not  on  its  dimension.  This  makes  a  systolic  design  ideal  for  implementing  finite 
difference  approximations,  where  the  number  of  mesh  points  is  typically  not'known 
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a  priori  but  where  the  maximum  bandwidth  is  determined  by  the  specific  scheme.  In 
the  next  section,  we  will  show  how  the  solution  of  (3.1)  can  be  accomplished  with 
systolic  arrays. 


CLOCK  PULSE  SYSTEM  CONFIGURATION 


2 


3 


4 


5 


6 


7 


Figure  3.2.4  A  walk  through  several  clock  cycles  of  the  systolic  layout. 
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3.3  Results  from  Linear  Algebra 

To  solve  (3.1),  we  use  Gaussian  elimination  without  pivoting.  This  method,  which  is 
stable  since  (/  —  At  An)  is  strictly  diagonally  dominant,  results  in  matrices  L  and  U 
such  that 

LU  ={I-  AtAn) 

where  U  is  an  upper  triangular  matrix  and  L  is  a  unit  lower  triangular  matrix. 

The  matrices  L  and  U  will  be  bi-diagonal  since  Gaussian  elimination  without 
pivoting  is  used.  As  is  standard  with  Gaussian  elimination,  once  the  factors  L  and  U 
are  found,  (3.1)  is  solved  by  finding  x  such  that 

Lx  =  DkV’‘  (3.3) 


and  then  by  solving 


(3.4) 


to  get 

Let  L  be  the  lower  triangular  matrix  resulting  from  the  factorization  of  a  band 
matrix.  We  are  interested  in  solving  Lx  =  b  for  the  vector  x.  That  is,  we  want  to 
solve 


11  0  \ 

/Xi  \ 

21  I22  0 

X2 

62 

31  ^32  hs  0 

X3 

— 

h 

0  I42  I43  liA 

X4 

b4 

'■.) 

[ :  J 

V :  y 

(3.5) 


This  is  solved  by  a  simple  back 


substitution  algorithm.  It  is  well  known  that  this  can 


be  accomplished  by  systolic  arrays  (Kung  &  Leiserson  [1980]). 

To  see  this,  notice  that  the  elements  Xi  can  be  computed  by  the  following  recur¬ 


sion: 


=  0, 

+  hkXk, 

Xi  =  {bi  - 


(3.6) 
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shows  the  overall  layout. 


Figure  3.3.2.  Overall  layout  for  the  sequential  detector. 
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3.4  Implementation  Issues 

In  collaboration  with  D.  Simmons,  a  member  of  the  research  staff  of  the  Systems 
Research  Center,  we  have  designed  a  board  to  implement  the  systolic  architecture 
discussed  above.  This  implementation  will  be  completed  shortly  and  will  reside  as  a 
special  purpose  processor  in  an  IBM  personal  computer.  The  hardware  details  of  the 
actual  design  and  layout  will  appear  in  a  technical  report.  In  order  to  get  some  idea 
of  the  implementation  constraints,  we  will  discuss  some  of  the  issues  raised  in  going 
from  the  block  diagram  layout  of  Figure  3.3.2  to  the  hardware  design. 

The  major  implementation  question  was  ‘  Should  we  use  fixed  point  or  floating 
point  arithmetic?’.  Among  the  advantages  of  fixed  point  arithmetic  are  that  fixed 
point  computations  are  up  to  four  times  faster  than  the  corresponding  floating  point 
and  that  fixed  point  arithmetic  leads  to  a  simpler  design.  However,  the  numbers 
in  our  problems  typically  have  large  dynamic  ranges.  It  is  common  to  have  density 
values  as  small  as  1  x  10~®  with  likelihood  values  as  high  as  1  x  10®.  In  addition, 
we  do  not  want  to  lock  ourselves  into  a  fixed  dynamic  range  as  would  happen  if 
we  implemented  fixed  point.  Therefore,  in  order  to  maintain  flexibility,  we  decided 
to  sacrifice  processor  speed  and  complexity  and  implement  the  design  using  IEEE 
standard  floating  point.  A  brief  descriptioon  of  the  processor  layout,  the  architecture 
and  the  components  of  proposed  board  prototype  design  are  shown  in  the  next  few 
pages.  We  have  named  the  proposed  chip  the  “Zakai  I”  chip.  Once  the  processor  is 
built  we  will  be  in  a  better  position  to  evaluate  the  optimal  design. 

Initial  timing  calculations  indicate  a  processor  speed  of  22  million  floating  point 
operations  per  second.  This  means  that  with  100  grid  points,  we  can  calculate  5000 
solutions  per  second. 


,  sets  1  and 


Zakai  Solver  Overview 


Zakai  1  Solver  Architecture 
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Actual  Systolic  Layout 


Inner  Product  &  Special  Processor 
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3.5  Summary 

We  have  shown  that  a  systolic  processor  can  be  used  to  implement  a  sequential  hy- 
popaper  tester.  It  was  seen  that  the  control  structure  for  the  resulting  processor 
was  simple  and  regular.  This  makes  it  ideal  for  VLSI  implementation.  Furthermore, 
timing  results  were  given,  hence  exact  computation  time  for  the  processor  can  be 
calculated. 

This  solves  the  sequential  hypothesis  testing  problem  for  the  scalar  x  case  unfor¬ 
tunately  more  work  remains  for  the  vector  x  case. 


4.  The  DELPHI  System 


We  are  currently  developing  an  expert  system  to  facilitate  the  CAD  of  sophis¬ 
ticated  and  complex  chips  for  non-linear  signal  processing.  The  software  system  is 
named  DEsign  Laboratory  for  Processing  Hidden  Information  (DELPHI)  and  can  be 
implemented  on  any  AI  machine  carrying  common  LISP. 

The  system  being  developed  combines  an  AI  engine,  symbolic  algebra  and  mul¬ 
tiprocessor  numerical  schemes.  It  has  the  capability  of  reasoning  mathematically  and 
makes  available  to  a  control  engineer,  sophisticated  mathematical  and  computational 
tools,  in  a  form  suitable  for  immediate  use.  Its  major  advantage  is  its  capability 
to  interact  symbolically  with  the  user.  The  block  diagram  of  the  DELPHI  software 
system  is  shown  in  figure  4.1  below.  The  system  can  be  used  as:  (a)  a  tool  for 
integrated  design,  (b)  an  advanced  teaching  aid  (we  have  successfully  used  it  with 
third  year  undergraduate  students),  (c)  a  tool  for  integrating  symbolic  and  numerical 
computation. 

The  “smart”  interface  block,  allows  the  user  to  enter  a  signal  and  observation 
model  symbolically.  We  are  currently  implementing  a  module  that  will  diagnose  user 
“expertness”,  so  as  to  allow  various  degrees  of  flexibility  to  the  user.  We  are  also 
implementing  a  module  capable  of  understanding  the  “nature”  of  the  module,  i.e., 
diffusion,  point  or  Markov  chain  type  models. 

The  modelling  block  can  currently  build  automatically  a  linear  dynamical  system 
model  for  a  Guassian  process  using  the  Akaike-Rissanen- Hannan  theory.  We  are 
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Figure  4.1  Blocks  of  the  DELPHI  System 
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currently  working  on  building  general  diffusion  type,  point  process  type  and  hidden 
Markov  chain  type  models.  This  block  will  have  the  ability  to  call  a  variety  of  statistics 
and  time  series  libraries  from  the  AI  machine  for  performing  statistical  validation  tests 
on  the  model  under  construction. 

The  likelihood  ratio  (or  sufficient  statistic  computer)  block  is  under  continuous 
development,  and  is  currently  the  most  advanced.  In  addition  to  the  results  of  the 
present  paper,  it  can  perform  similar  computations  for  point  process  observations, 
mixed  diffusion  point  process  type  models,  and  Markov  chain  models. 

During  next  year  we  will  couple  the  system  to  a  silicon  compiler  for  actual  VLSI 
chip  layout.  What  we  will  really  develop  is  a  “smart”  compiler  which  can  preselect 
the  architecture  to  best  fit  the  problem  based  on  the  parallel  complexity  theory  of  the 
Zakai  equation. 

The  likelihood  ratio  block  has  currently  the  following  capabilities:  (an  M  indi¬ 
cates  a  MACSYMA  computation,  an  F  indicates  a  FORTRAN  computation). 

Ml:  Input  f,g,h  in  symbolic  form.  Generate  the  multidimensional  Zakai  equation. 
M2:  Compute  automatically  discretizations  of  all  stochastic  differential  equations. 
Automatically  generate  FORTRAN  code. 

F3:  Solve  numerically  pathwise  the  resulting  stochastic  difference  equations  and  store 
the  y-  paths. 

M4:  Generate  discretization  schemes  for  the  Zakai  equation  and  automatically  gener¬ 
ate  the  associated  FORTRAN  code. 

F5:  Automatically  integrate  numerically  the  discretized  Zakai  equation. 

F6:  Display  graphically  the  solution. 

M7:  Compute  symbolically  discretizations  of  the  Likelihood  Ratio  and  automatically 
generate  appropriate  FORTRAN  code. 

F8:  Evaluate  numerically  the  likelihood  ratio. 

Further  additions  currently  planned  include:  apriori  bounds  for  performance  of 
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detectors  and  estimators,  additional  numerical  schemes  for  the  Zakai  equation  (e.g., 
multi-grid  algorithms),  automatic  performance  evaluation  of  sequential  detectors.  In 
addition  we  have  initiated  the  development  of  a  LISP  based,  domain-specific,  higher 
level  language  for  signal  processing.  Finally,  a  reduced  version  of  the  system  is  being 
ported  on  an  IBM  PC  AT. 


Future  Directions 


We  have  presented  a  VLSI  design  to  implement  a  real  time  processor  for  the  case 
where  the  signal  process  is  scalar  valued.  More  research  work  needs  to  be  done  in 
order  to  discover  the  appropriate  architecture  for  the  general  case  where  the  signal 
process  is  multidimensional.  It  is  clear  that  this  will  be  accomplished  by  exploiting 
problem  specific  information  into  the  design.  It  is  also  clear  that  an  implementable 
architecure  will  be  very  regular  but  exceedingly  redundent  and  detailed.  Therefore, 
a  systematic  and  automatic  method  must  be  devised  which  takes  .a  problem  and 
produces  a  VLSI  design. 

It  appears  that  systolic  arrays  can  be  used  for  higher  dimensions,  up  to  about 
5  or  6.  For  even  higher  dimensions  it  is  clear  that  we  will  need  massively  parallel 
architectures. 
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